Optical excitations in organic molecules, clusters and defects studied by 
first-principles Green's function methods 
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Spectroscopic and optical properties of nanosystems and point defects are discussed within the 
framework of Green's function methods. We use an approach based on evaluating the self-energy 
in the so-called GW approximation and solving the Bethe-Salpeter equation in the space of single- 
particle transitions. Plasmon-pole models or numerical energy integration, which have been used in 
most of the previous GW calculations, are not used. Fourier transforms of the dielectric function 
are also avoided. This approach is applied to benzene, naphthalene, passivated silicon clusters 
(containing more than one hundred atoms), and the F center in LiCl. In the latter, excitonic effects 
and the Is — > 2p defect line are identified in the energy-resolved dielectric function. We also compare 
optical spectra obtained by solving the Bethe-Salpeter equation and by using time-dependent density 
functional theory in the local, adiabatic approximation. From this comparison, we conclude that 
both methods give similar predictions for optical excitations in benzene and naphthalene, but they 
differ in the spectra of small silicon clusters. As cluster size increases, both methods predict very 
low cross section for photoabsorption in the optical and near ultra-violet ranges. For the larger 
clusters, the computed cross section shows a slow increase as function of photon frequency. Ionization 
potentials and electron affinities of molecules and clusters are also calculated. 



I. INTRODUCTION 

Universality and predictive power are characteristic 
features of good first-principles theories. Such theories 
are invaluable as guides for experimental work, espe- 
cially when knowledge about the system being inves- 
tigated is limited. Density-functional theory (DFT) is 
recognized as state-of-the-art theory for the calculation 
of various ground-state properties of electronic systems 
from first principles, having been applied to a variety 
of different systems ranging from crystalline semicon- 
ductors to surfaces and nanostructuresi In contrast, 
DFT has known difficulties in predicting quantities asso- 
ciated with excited states, which is not surprising since 
it was first proposed as a ground-state theory. As an ex- 
ample, the electronic band structure provided by DFT 
consistently shows an underestimated gap between va- 
lence and conduction bandsi 2 *^. As consequence, the 
electronic band structure provided by DFT gives poor 
predictions for the onset of photoemission and inverse 
photoemission 2 ^. Exact-exchange functionals improve 
results considerably for molecular systems^, but more 
work is needed in crystalline systems^. Gap narrowing 
also affects the linear optical spectrum, when calculated 
from single-electron transitions from the valence band to 
the conduction band^i 

Spectroscopic and optical properties of semiconductors 
have been successfully predicted from ab initio Green's 
function methods^ 7 -. The GW approximation has been 
used to provide the quasiparticle band structure of semi- 
conductors, large-gap insulators, metals, surfaces, nanos- 
tructures and other materials^. By solving the Bethe- 
Salpeter equation for electrons and holes, the linear op- 
tical spectrum can also be obtained, providing reliable 



results for the optical gap and excitonic effectsi&SiiSiii. 

Since knowledge of the electron self-energy is required in 
solving the Bethe-Salpeter equation (BSE), this method 
is often refered to as GW+BSE^. Although powerful, 
the GW+BSE method is very complex when applied to 
non-periodic systems, especially when discrete Fourier 
analysis is not applicable and approximations such as 
the generalized plasmon pole models do not lead to sub- 
stantial simplification. Some materials that fall in this 
category are confined systems: quantum dots, clusters, 
molecules, nanostructures. On the other hand, nanos- 
tructures are now the subject of intense work due to 
promising technological applications and recent advances 
in the preparation of artificial nanostructuresi&ii At the 
same time, significant effort has been dedicated on the 
theory side. Some of the difficulties one may encounter 
by applying the GW approximation to large molecules 
were emphasized recentlj*i^ii&. So far, reliable calcula- 
tions of electronic self-energy within the GW approxima- 
tion have been done for confined systems with up to a few 
tens of atoms^iiLiSiiS. Similarly, the Bethe-Salpeter 
equation has been solved for a limited number of con- 
fined systems. During the last decade, optical properties 
of confined systems have been investigated within time- 
dependent DFT, often with very good result o 7 i 20 . Time- 
dependent density DFT in the local, adiabatic approxi- 
mation (TDLDA) is formally simpler than the GW+BSE 
theory, but at the same time directed towards solving 
the same problem: predicting optical properties and neu- 
tral excitations in an electronic system. A drawback of 
TDLDA is that it reportedly fails to predict correct opti- 
cal gaps and excitonic effects in extended system o 7 ! 21 . Al- 
ternative TDDFT functionals which include many-body 
effects have been proposed 21 i 22 i 23 i 24 i 25 . 
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As mathematical simplification, applications of the 
GW+BSE method in nanostructures frequently impose 
artificial periodicity by placing the electronic system in a 
large supercelUt This is a perfectly justifiable procedure, 
but it has the obvious deficiency of increasing the com- 
putational effort when supercells are required to be large. 
A more efficient procedure would be to take advantage of 
properties intrinsic to confined systems and design a nu- 
merical implementation of such Green's function meth- 
ods, and we now propose one such implementation. A 
key ingredient is electronic screening. We incorporate it 
into the theory by computing the electron polarizability 
function in two approximations: random-phase approxi- 
mation (RPA) and TDLDA. The electron self-energy is 
calculated from first principles and the Bethe-Salpeter 
equation is set and solved in the space of single-particle 
transitions. In this framework, many-body functions are 
expressed as matrices in this transition space and in en- 
ergy (frequency) representation. Fourier transforms are 
not needed, and energy integrations are evaluated by in- 
tegrations over poles. Portions of this work have been 
reported earlier—. 

The methodology proposed here is applied to a number 
of interesting cases. One of them is isolated oligoacenes, 
for which extensive experimental data is available but 
the only GW-based analysis known to us have been pub- 
lished recently^. The second application is in silicon 
nanoclusters, studied previously within TDLDA2S and, 
for the smaller clusters, GW+BSE£ii&. Here, we investi- 
gate clusters with bulk-like crystalline structure and re- 
port benchmark GW+BSE calculations for clusters with 
more than one hundred silicon atoms. The last applica- 
tion is in the F-center defect in LiCl. This is a challenging 
case, with characteristics of both confined and periodic 
system. We show that the resulting optical spectrum cor- 
rectly contains both the signature of defect levels located 
within the band gap and electron- hole interactions. The 
latter manifest themselves in the energy-resolved dielec- 
tric function as an exciton feature below the electronic 
band gap. TDLDA is shown to predict an optical tran- 
sition between two defect levels, but the exciton feature 
is not found. The paper is organized as follows: in Sec- 
tion [HJ we present the theoretical framework, starting 
from TDLDA and continuing through the GW method 
and solution of the Bethe-Salpeter equation. Section ILTTI 
is devoted to applications, containing one subsection for 
each particular system: oligoacenes, silicon clusters, and 
F center in LiCl. We draw final conclusions in Section 

El 



II. THEORY 
A. Linear response within TDLDA 

Gross and KohnSi have shown how to extend DFT 
to the time-dependent case by analyzing the effect 
on the charge density upon the action of an exter- 
nal potential that changes in time. Their theoret- 
ical approach can be further simplified by making 
use of two assumptionsii22i2&: adiabatic limit, and 
local-density approximation, under which the exchange- 
correlation kernel is instantaneous in time and a sole 
function of charge density at each point in space. 
This is the time-dependent local-density approximation 
(TDLDA) 20 ! 28 ! 29 ! 30 . 

When the external potential is due to an applied elec- 
tric field, the linear response in charge density p is related 
to the external potential via a response function (the po- 
larizability) according to 



n/(i,2) 



SV ext {2) 



(1) 



where we use a many-body notation for space-time and 
spin variables: (1) = (ri, tx,rx). The subscript / in- 
dicates that the response function above is calculated 
within TDLDA. Working in frequency domain^, the 
TDLDA polarizability can be written as a sum over nor- 
mal modes, 



n/(r,r';£) 



2E s ftWft(r') 
l l 
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(2) 



where the normal modes of excitation in the sys- 
tem are denoted by uj s , and + represents a positive 
infinitesimal^ 3 . . We assume Ti = 1 throughout. The factor 
of 2 comes from summation over spin indices. Still fol- 
lowing the frequency-domain formulation by Casida™*^, 
the amplitudes p s (r) are obtained by solving a gener- 
alized eigenvalue problem. In order to obtain the de- 
sired eigenvalue equation, we first expand p s in a series 
of single-particle transitions from an occupied level v to 
an unoccupied level c: 



Ps{r) = ^X s vc y v {r)ip c (r) 



1/2 



(3) 



where Kohn-Sham eigenvalues are denoted Si, with cor- 
responding cigenfunctions ip{. For clarity, we reserve in- 
dices v,v' for occupied orbitals and indices c,c' for un- 
occupied orbitals. The coefficients X above satisfy an 
eigenvalue equation in (uj 2 )2^: 



R 1 / 2 [R + 4(K* + K LDA )] R^ 2 X = w s 2 X , (4) 
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where R, K 31 , and K. LDA are matrices in the space of 
single-particle transitions: 



(5) 



K% cv , e = / dr / drV,(r)^ c (r)y(r,r')^(r')^(r') 



(0) 



K$% A d = I dr^(r)^ c (r)/ xc (r)9v(r)¥v(r) • (7) 



The generalized eigenvectors X are normalized so that 
X S X S = 5 SS >. Throughout this article, eigenfunctions 
<p are assumed to be real functions. Indeed, they can 
be made real if the Kohn-Sham Hamiltonian is real 
and the electronic system is confined. This assump- 
tion is used in Eqs. |J2}, © and (@J. A notable situ- 
ation when the equations above should be modified is 
in periodic systems, when eigenfunctions are often ex- 
pressed as Bloch functions. The generalization to com- 
plex eigenfunctions is knowniSi. Also, we assume that 
the system has an energy gap and it is spin-unpolarized, 
although this is not an essential assumption and the 
same framework holds for gap-less systems (e.g, the F 
center in LiCl, Section IIII C|) as well. There are ex- 
tensions of TDLDA to spin-polarized and/or gap-less 
systems2S^2*^i. Exchange-correlation effects are included 
in the kernel f xc = which is a local, energy- 

independent quantity within TDLDA. The Coulomb po- 
tential is simply V(r,r') = rgzyi ■ 

Once the eigenvalue problem is solved, one can easily 
compute the electrostatic polarizability tensor, 



a 



(11) 



Fundamental quantities in linear-response theory are 
the (longitudinal) dielectric function e and its inverse 
e _1 A2i. In order to discuss them, we remind ourselves of 
the Kohn-Sham system: consider a set of non-interacting 
electrons subject to a Hartree potential and an exchange- 
correlation potential V xc , chosen so that the charge den- 
sity of this fictitious system and of the real system are 
identicali*S. The presence of an external potential in- 
duces charge redistribution and polarization in the Kohn- 
Sham system, which results in partial screening of the ap- 
plied potential SV ex t- For a time- varying external pertur- 
bation, electrons are subject to an effective perturbation 
potential given by 

SV eff [p](l) = SV ext (l) + 6Vscf[p}(1) , 

where the self-consistent field is affected indirectly via 
charge density, SVscf = VSp + ^^Sp^. Defining the 
inverse dielectric function as the change in effective po- 
tential due to an external perturbation^^, 



e- 1 (l,2) = 



SV ext (2) 



we then obtain a relationship between inverse dielectric 
function and polarizability: 



e~ 1 (l,2) = 6(1,2)+ /(3)d(3) 



V(l,3) 



6V XC (1) 



6 P (3) 



11(3,2) 



In frequency representation, and using TDLDA for the 
exchange-correlation potential V xc , the above equation 
reduces to 



e 2 Tp 



(8) 



where m is the electron mass and J-@ 7 is the oscillator 
strength along the three cartesian components /3, 7 = 
{x,y,z}, 



= 4mu>. (y dvp s (r)(3^j (J drp s (r) 7 



(9) 



By construction, the TDLDA polarizability satisfies the 
oscillator-strength sum rule, 



\^ = (number of valence electrons) 



(10) 



Another quantity of interest is the cross section for 
absorption of light in a confined system: 



ej\r, r'; E) = 5(r, r') + / dr" [V(r, r") 

+ f xc (r)6(v-r")}IL f (r",r';E) . (12) 

Along similar lines, the dielectric function e is given by 

e/(r, r'-E) = S(r, r') - / dr" [V(r, r") 

+f xc (v)5(v-v")]xo(r",r';E) , (13) 

where \o is the irreducible polarizability operator, within 
the random-phase approximation (RPA)2LS: 



Xo(r, r'; E) = 2J2 S ^(r)<^ c (r)^,(r')(^ c (r') 



E-E c +Ev+iO+ E+E c -Ev-iO + 



(14) 



Eqs. I|12fl an( l l|13|) describe screening of the exter- 
nal scalar field and of the induced field. Recently, Eq. 
(I13fl was used to calculate the TDLDA static screening 
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in silicon clusters from first principles*^. The importance 
of Eq. H12fl is that it provides a direct relationship be- 
tween polarizability and inverse dielectric function. For 
completeness, we present expressions for the dielectric 
function and its inverse within the RPA: 



e - 1 (r,r'; J B)=5(r,r')+ / dr"V(r, r")Ho(r", r'; E) 



(15) 



X (l,2)w-iG(l,2 + )G(2,l)= Xo (l,2) . 

The screened Coulomb interaction is evaluated in terms 
of the dielectric function. From Eq. I|17[l . one gets: 

e (l,2)=(5(l,2)- /d(3)y(l,3) X o(3,2) , 



e {v,v';E) = 5{v,v')- j dr'V(r, r")xo(r", r'; E) , 

(16) 

where IIo is the RPA polarizability, evaluated from Eqs. 
©, © and (QJ after setting f xc = 0. 



Wo(l,2) 



Jd(3)e - 1 (l,3)n3,2) 



= V(l,2) + J d(34)y(l,3)n (3,4)y(4,2) 
and the self-energy is 



B. Electronic self-energy 

A key quantity in solving the Bethe-Salpeter equa- 
tion for optical excitations is the electronic self-energy 
S, which can be computed from first principles within 
the GW approximation (GWA) 2 i 3 i 4 . The basis of this 
approximation lies in the so-called "Hedin equations", 
a set of non-linear many-body equations which relate 
the self-energy with Green's function G, polarizability 
X, screened Coulomb potential W, and vertex function 

14 



E(l,2) = iG(l,2)W (2,l+) 



(21) 



W(1,2) = V(1,2)+ / d(34)y(l,3)x(3,4)?y(4,2) 



(17) 



2) = -i / d(34)G(l, 3)G(4, l+)r(3, 4; 2) , (18) 



(19) 



E(l,2) = i J d(34)G(l,3)VF(4,l+)r(3,2;4) , 



r(l,2;3) =5(1,2)5(1,3) 

+ / d(4567) Jggg G(4, 6)G(7, 5)r(6, 7; 3) (20) 

The approach taken by Hedin essentially generates a 
perturbation series in the screened interaction W. This 
expansion is expected to converge faster than an ex- 
pansion in powers of the bare Coulomb potential V, as 
long as electronic screening is strong. Following this 
assumption**, the self-energy in Eq. (12011 is first taken 
as zero and the vertex function reduces then to a delta 
function: 

r(l, 2; 3)«i(l, 2)6{1, 3) . 

With that vertex function, the polarizability x reduces 
to the RPA**: 



The first ab initio calculations of self-energies for 
semiconductors*****^ followed the method outlined above. 
For periodic systems, the dielectric function can be con- 
veniently expanded in plane wave basis and numerically 
inverted. Matrix inversion is a reasonably inexpensive 
step when the plane- wave expansion has a few hundreds 
or thousands of components. Dynamical effects can be 
included for instance by inverting the dielectric function 
at some values of frequency and either performing nu- 
merical integration over the frequency axis«*>i***, or using 
a generalized plasmon pole model**. 

On the other hand, numerical inversion of the dielectric 
matrix in real space is problematic because of the size of 
the matrix. As an example, a converged calculation for 
the ground state of the silane molecule (S1H4) requires a 
real-space grid containing about 10 5 points****. Inverting 
the dielectric function in real space would require invert- 
ing a dense matrix of size 10 5 x 10 5 . Although straight- 
forward, this task requires an extremely large amount of 
floating point operations and CPU memory. This is cer- 
tainly an extreme situation, but it illustrates the type of 
problems involved with straight matrix inversion. The 
same issue is present in periodic systems with large pe- 
riodic cells: since the number of plane-wave components 
in a Fourier expansion is proportional to the volume of 
the cell, direct matrix inversion of e is also numerically 
expensive. 

Significant numerical simplification is achieved 
by working in the representation of single-electron 
transitions****, instead of using plane-wave or real-space 
representations for the dielectric function, there are 
two major advantages in doing that: the space of 
transitions is often much smaller than either the real 
space or reciprocal space, leading to matrices of reduced 
size; and integrations over frequency can be performed 
analytically since the polarizabilitics xo and Ho have 
known pole structure. In the space of transitions, a 
matrix element of the self-energy between Kohn-Sham 
orbitals j and j' is given schematically by: 
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(mE')\j') = J ^ e ~ lE0+ G(E'-E) [V + Vn Q (E)V] 

(22) 

The integral above can be replaced with a sum over poles 
below the real energy axis. Following Hcdin^i, we write 
it as a summation of two terms: a bare exchange contri- 
bution H x and a correlation contribution £ c , 

occ. 

(j|£*|j'} = -£^W , (23) 



m c (E)\f) = E-7-lr, ' (24) 

71 S 

where 

^■=£^c(^y /2 ^ c . (25) 

vc ^ ' 

The derivation of Eqs. H3J), (JUJ) and JUJ) from Eq. 
(|22|l is outlined in Appendix ^] Whereas the summa- 
tion over n in Eq. (|23|l is performed over all occupied 
Kohn-Sham orbitals, Eq. (|24|l has a summation over 
occupied and unoccupied orbitals. One can accelerate 
convergence in the last summation by truncating it and 
evaluating the remainder under the static approximation 
(see Appendix |B|) . The coefficient r/ n has value +1 for 
unoccupied orbitals and -1 for occupied orbitals. In the 
equations above, eigenvectors X and eigenvalues to are 
still evaluated under the random-phase approximation, 
obtained by setting f xc = in Eq. In the following 
discussion, we refer to this level of approximation for the 
self-energy as GWo- 

At this point, it is natural to advance one step fur- 
ther and use dielectric screening within TDLDA instead 
of RPA. The impact of TDLDA screening in the con- 
text of the GW approximation in bulk silicon has been 
analyzed before^iSiiS. In order to have a controlled 
level of approximation in Hedin's equations, we return 
to Eq. (|2t)(l and start an iterative solution by assum- 
ing £(1,2) w V XC {\)5{1,2). The vertex function then 
becomes 

r(l,2;3) ra 6(1, 2)8{2, 3) + / d(45) x 

[-iS(l, 2)/ xc (l)] G(l, 4)G(5, l+)r(4, 5; 3) (26) 



The irreducible polarizability \ is now 



X (l, 2) = X o(l, 2) + J d(3)xo(l, 3)/xcX(3, 2) . (27) 

For the screened Coulomb interaction, we write Eq. 
(|17(l in terms of the full polarizability operator: 



W(1,2)=V(1,2) + J d(34)V(l, 3)11(3, 4)V(4, 2) , 
with 

n(l,2)= X (l,2)+ /d(34) X (l,3)V(3, 4)11(4, 2) . 
From Eq. ifSTjl. we get 

n(l, 2) = X o(l, 2) + / d(34)xo(l, 3) [F(3, 4) 
+ f xc (3)S{3, 4)] IT(4, 2) . 

Although written in real-space representation, the 
function II above is identical to Eq. J5J), and we refer 
to it as II f in the subsequent discussion. This function 
describes polarization due to an external potential ac- 
companied by dynamical screening produced by the self- 
consistent field. Finally, we arrive at the following ex- 
pression for the self-energy: 

£(1,2) =iG(l,2) [V(l+,2)+ fd(34){y(l,3) 

+ / ac (l)5(l,3)}n / (3,4)F(4,2)] . (28) 

We note that the use of TDLDA screening in the self- 
energy causes the inclusion of a vertex term (the second 
term inside curly brackets above) . This level of approx- 
imation has been used before in the study of the quasi- 
particle band structure of crystalline silicon^. Eq. I|28() 
is a valid approximation for the self-energy but it is not 
symmetric with respect to the interchange of indices 1 
and 2. In principle, this symmetry is not present in Eq. 
(|H1 . but it can be recovered by defining a "left-sided" 
vertex function written schematically as £ = iTWG as 
opposed to £ = iGWT (c.f. Eq. 11^ . This deficiency 
is corrected by symmetrizing the last term in Eq. 128|l . 
Schematically, we rewrite £ = iG[V + VU f V + fUfV] as 
Y t = iG[V + VUfV+lVn f f+yUfV]. Matrix elements 
of the self-energy are then given by: 



(mE')\f) = w 



dE 
< 

2tt 



G(E' - E) 



V 



VU f (E)V + -VU f (E)f xc 



;f,cU f (E)V 



\f) 



(29) 



This self-energy operator is written as a sum of three contributions: £ = ~E X + £ c 



£y, where the first two 



6 



(a) 




FIG. 1: Feynman diagrams for the polarizability operator IT 
(a,b) and the function $ (c,d). The local exchange-correlation 
kernel is represented by a crossed line. Solid oriented lines 
are Green's functions. Dashed lines denote the bare Coulomb 
potential. 

terms arc given by Eqs. I|23|) and i|24[l respectively, but 
now with full TDLDA screening ( f xc ^ 0). The last term 
is a vertex correction: 

, , , , V s F s -, + F s -V s -, 

(iis / (^)i/) = EE ZZ _7 " 3 . ( 3 °) 

with 

/ — \ 1 / 2 

ps _ K LDA I £c e _v_ \ /o-i^i 

n J / j njvc I ^ J vc ' \ / 

vc V * / 

In order to make distinction between the two levels of 
approximation, we refer to the last approximation (Eq. 
E| as GW/, as opposed to GW , Eq. (J22J)- 

It is not unexpected that, by using a polarizability 
function from TDLDA, we obtain a self-energy operator 
that has vertex corrections. In the language of many- 
body physics, the assumption in Eq. i|26|> implies that 
additional Feynman diagrams are included in the po- 
larizability and in the self-energy. Fig. ^a,b) shows 
the diagrams included in the polarizability II for both 
RPA and TDLDA. Although such diagrams include the 
LDA kernel, which does not have expansion in G and 
V, we can still define a many-body function $ such that 
£(1,2) = 6$/<5G(l,2). The diagrams for $ are depicted 
in Fig. ^c,d). Adding the last diagram in Fig. ^b) 
amounts to enhanced screening, since it allows for a new 
channel for electrons to redistribute in the presence of 
an external potential. This is counterbalanced by the in- 
clusion of a vertex diagram (right-most diagram in Fig. 
^d)). We also note that Fig. provides a justification 
for symmetrizing Eq. I|28[) : differentiating the additional 
vertex diagram with respect to G leads to a pair of dia- 
grams consistent with Eq. (|3U|) . 



Inclusion of self-energy corrections in the description 
of the electronic system is done following the usual pro- 
cedure: we assume the quasiparticle approximation and 
write down an eigenvalue equation for quasiparticles^, 

[H LDA + £ - V xc \ V>,- = Eji/fj , (32) 

where quasiparticle orbitals ipj are expanded in the ba- 
sis of Kohn-Sham cigenfunctions and Hlda is the (di- 
agonal) LDA Hamiltoniani. Quasiparticle energies and 
wave-functions are now found by direct diagonalization 
of the equation above. 

In writing Eq. I|32|) . we should consider carefully the 
energy dependence of £. Hybertsen and Louia^ have 
shown that evaluating the self-energy around the quasi- 
particle energy leads to accurate band structures for cu- 
bic semiconductors. Since the quasiparticle energy is 
not known before Eq. (|32[) is solved, the suggested 
procedure^ is to evaluate the operator (T,(E) — V xc ) and 
its energy derivative at the LDA eigenvalue, E = Eld a, 
and use linear extrapolation for the actual quasiparticle 
energy. We follow a similar methodology: first, we evalu- 
ate the diagonal part of the quasiparticle Hamiltonian in 
Eq. (|32f> and obtain a first estimate for quasiparticle en- 
ergies by solving the equation Ej = Ej + (j\Y,(Ej) — V xc \j). 
In the next step, we include off-diagonal matrix elements 
and proceed through full diagonalization. At the end, 
quasiparticle energies are still close to their first estimate. 

An open issue regarding the GW method is whether 
self-consistency between self-energy, polarizability and 
Green's function should be imposed or not and, if so, 
how to do i1w When the GWo and GW/ approxi- 
mations were obtained by iterating Hedin's equations, 
self-consistency was lost and the vertex function was 
drastically simplified. Early work has indicated that 
self-consistency and vertex contributions partially can- 
cel each othesafi, which may explain the remarkable suc- 
cess of the GW method for semiconductors. There have 
been some attempts at imposing partial self-consistency 
between self-energy and Green's function^*, which re- 
sulted in quasiparticle energies improved by a fraction 
of electron-volt. Nevertheless, full self-consistency be- 
tween self-energy and Green's function was observed to 
degrade the quasiparticle bandwidth and the description 
of the satellite structure in the electron gas^i. In the sub- 
sequent applications, we do not attempt to impose self- 
consistency. Instead, the single-particle Green's function 
is always constructed from Kohn-Sham eigenvalues and 
cigenfunctions. 



C. Bethe-Salpeter equation 

While the GW method provides a description of quasi- 
particles in the system (i.e., the energy needed to add or 
extract one electron from the system), the solution of the 
Bethe-Salpeter equation gives information about neutral 
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excitations: the process of promoting electrons from oc- 
cupied quasiparticle orbitals to unoccupied ones. The 
important Green's function now is no longer the one- 
electron Green's function but the two-particle or, more 
specifically, the electron-hole Green's functional. 

We start with a many-body expression for the polariz- 
ability: 



11(1,2) = -iL(l,2;l+2 H 



(33) 



where L is the electron-hole correlation functionLSiSi, 
which satisfies the Bethe-Salpeter equation: 

L(1,2;3,4) = G(1,4)G(2,3) 
+ / d(5678)G(l,5)G(6,3)A'(5,7;6,8)i(8,2;7,4) (34) 

The kernel operator K describes interactions between 
the excited electron and the "hole" left behind in the 
electron sea. The connection with GW is present in two 
aspects: 

1. The Green's function that enters in Eq. H34|) is 
the Green's function for the interacting electronic 
system. Within the quasiparticle approximation, it 
is calculated with eigenvalues and wave-functions 
obtained by diagonalizing Eq. I|32|l . 

2. The kernel K is related to the self-energy byi2*2i: 



tf(l,2;3,4) = -W(l,3)<5(2,4)V(l,2) + 



<5£(1,3) 



(35) 



Although the functions G and K can be evaluated in 
many approximations, it is important to retain the same 
level of approximation in both quantities, so that all 
important Feynman diagrams are included and no dia- 
grams are double-counted. In the lowest level of approx- 
imation, many-body effects can be ignored completely 
and the self-energy approximated by the LDA exchange- 
correlation potential, £(1,2) w V xc (l)8(l, 2). In this 
case, the Green's function to be used is constructed from 
LDA eigenvalues and eigenfunctions, and Eq. (|34|) re- 
duces itself to the TDLDA eigenvalue equation, Eq. 10} . 
This fact has been explored in recent studies where model 
exchange-correlation functionals are designed to work as 
good approximations for the self-energy in Eq. I|35|l while 
retaining as much as possible the formal simplicity of 
TDLDAi2L2^. By its own nature, this family of approx- 
imations does not have Feynman diagrams in terms of 
which the self-energy is expanded. 

Another level of approximation is G Wo , Eq. I122II . The 
kernel K has then two terms: a bare exchange interac- 
tion, originated from the first term in Eq. I|35|) . and a 
screened interaction, from <5£/<5G = iWq. This approxi- 
mation has been used in a variety of different applications 
with very good results (see e.gm&). 



Finally, a third level of approximation is GW/, Eq. 
(|29|l . At this level, the kernel has an additional term 
besides the previous two: a LDA vertex correction. As 
opposed to the kernel present in the TDLDA equation, 
the BSE kernel within either GWo or GW/ has explicit 
energy dependence, from the polarizability operator n. 
This dependence was observed to be negligible when the 
interaction kernel itself is weakS. Ignoring dynamical ef- 
fects and ignoring the mixing of absorption and emission 
contributions in the kernel^, one can rewrite Eq. I|34|) 
as an eigenvalue equation^: 



n x A l vc = {E c -E v )A l vc 

+ E»'c' {^^vcv'c' + K tcv'c> + Klcv'S) A l v > c > , (36) 

where E Cl E v are quasiparticle energies of unoccupied and 
occupied orbitals respectively. The exchange kernel K x 
has the same form given by Eq. © . The kernels K d and 
are given by: 



V s , V s , 

K d ,, = K x , I \ > m cc 

vcv c' ^-^vv'cc' 1 / j 



ys ps , ps ys 
r vv' cc' 1 vv' cc' 

UJ S 



(37) 



(38) 



The eigenvectors A are normalized in the usual way, 
and the polarizability is now given by 



IL BS E(r,r';E)=J2 



Pijr)pi(r') _ P i{v)pi{t') 
E-Qi + iO+ E + fli~ i0+ 



with 



Pi 



(39) 



(40) 



As in the TDLDA, the solution of the Bethe-Salpeter 
equation provides the electrostatic susceptibility tensor 
and absorption cross section by equations analogous to 
(I5|).lj§|l. (|llfl . with the replacements p s — > pi and uj s — ► 0;. 



D. Technical considerations 

Working in the space of single-particle transitions has 
two clear advantages: it does not assume any particu- 
lar boundary condition over electron wave-functions, be- 
ing equally valid for confined as well as extended sys- 
tems; and it often leads to matrices smaller than the 
corresponding ones in real-space or reciprocal-space rep- 
resentations. In addition, the necessary numerical tasks 
are simple: matrix algebra (diagonalizations and matrix 
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products), and evaluation of integrals of the type in Eqs. 
© and Q . We can also take advantage of existing point 
symmetries in the system at hand and divide the transi- 
tion space into different group representations, resulting 
in block-diagonal matrices. This method can be readily 
extended to periodic systems by adding lattice periodic- 
ity as an additional symmetry operation. Each block is 
then associated to a particular k-point in the Brillouin 
zone. For large systems, the solution of Eqs. J3J and 
(|36[l may involve large matrices, but further simplifica- 
tion is obtained by restricting the transition space to low- 
energy transitions. In the calculation of the self-energy 
terms S c and £/, one can proceed through the sum over 
single-particle states by including the static remainder 
(see appendix El . 



III. APPLICATIONS 

In order to assess the reliability of the current imple- 
mentation, we have employed it to study the excited- 
state properties of different classes of systems, both lo- 
calized and extended. In all applications, we start from 
the electronic system in its ground state, as described 
within DFT. The exchange-correlation potential is used 
in the Ceperley-Alder formii With the use of norm- 
conserving pseudopotentials^, we are able to remove core 
electrons from the problem, reducing considerably the 
numerical effort. In all applications, we use pscudopo- 
tcntials constructed according to the Troullicr-Martins 
prescription^. We employ a real-space approach in the 
solution of the Kohn-Sham Equation^. In this approach, 
wave-functions are expressed directly as functions of po- 
sition. For a finite system, wave-functions are required 
to vanish outside a spherical boundary with adjustable 
radius. The volume inside this boundary is sampled by 
a homogeneous grid with fixed grid spacing h along the 
three Cartesian directions. Results are tested for conver- 
gence with respect to the grid spacing and the boundary 
radius. We explore point symmetries both in the solution 
of the Kohn-Sham equation and in the calculation of in- 
tegrals in Eqs. © and JJJ). For the self-energy operator, 
we include the static remainder according to appendix 151 
We carefully test convergence with respect to the high- 
est state explicitly included in the n summation. Ex- 
cept when noted, we compute self-energies at the GW/ 
level. Numerical accuracy in self-energy matrix elements 
is 0.2 eV or better. Excited states of the electronic sys- 
tem are calculated without the inclusion of effects due to 
structural relaxation. 



A. Benzene and naphthalene 

Having a relatively small number of valence electrons 
makes the benzene molecule, CqHq, a good test system 
for the theory just presented. In addition, benzene has 
been the subject of extensive experiments, probing vari- 



> 

S-H 

w 

N 



25 



20 



a is 

o 



v 
o 

-a 
H 



10 



- 


1,111 

C 6 H 6 


! | i i i i | i i i i | y 




❖ LDA 
X GW f 


/ o _ 


- 


+ GW o 


/* o 


- 




±3r o 
o 


_ 


w ❖ 
o 

i i i i i 






10 


15 20 25 



Experimental Ionization Energy (eV) 

FIG. 2: Ionization energies of benzene, associated to all oc- 
cupied orbitals in the molecule, as predicted within GWo and 
GW/. The theoretical ionization energy is the negative of a 
quasiparticle energy eigenvalue. Experimental data fror 



ous properties such as linear optical response, electronic 
structure and ionization^. On the theory side, its optical 
properties have been investigated with both TDLDA and 
GW-BSE approaches22i2£i2£. Quasiparticle energy levels 
of benzene and other oligoacenes have also been recently 
calculated using a first-principles Gaussian orbital-based 
method as well as a tight-binding approach^. The se- 
ries of oligoacenes, benzene being the first element, is a 
sequence of molecules composed of n co-joined aromatic 
rings, starting from n = 1 (benzene), n = 2 (naphtha- 
lene), etc. Oligoacene crystals have received renewed at- 
tention recently due to their potential use in organic thin 
film devicesii*^. 

In our calculations, we solve the Kohn-Sham equa- 
tions on a regular grid with grid spacing 0.4 a.u. (the 
Bohr radius being 1 a.u.). Electronic wave-functions 
were required to vanish outside a sphere of radius 16 a.u. 
centered on the molecule. Carbon-carbon and carbon- 
hydrogen bond lengths were fixed at their experimental 
values. The calculation of the self-energy operator was 
done including up to 305 unoccupied (virtual) orbitals in 
the TDLDA polarizability and in the calculation of cor- 
relation and vertex parts. This ensures that the sum rule 
is satisfied to within 30%, which was found to be suffi- 
cient for an accuracy of 0.1 eV in quasiparticle energies. 
After the self-energy operator was computed, quasipar- 
ticle energies for occupied and unoccupied orbitals were 
obtained by diagonalizing the Hamiltonian in Eq. I|32|) . 
Extensive convergence tests have been performed on all 
relevant numerical parameters: choices of 0.3 and 0.4 a.u. 
were used for the grid spacing; we have used boundary 
radii of up to 20 a.u., and included up to 1000 virtual 
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orbitals. For the last choice, the sum rule is satisfied to 
within 14%. 

Fig. |21 shows a comparison between theoretical and 
experimental ionization energies. In the framework of 
the GW approximation, the first ionization energy is in- 
terpreted as the negative of the quasiparticle energy for 
the highest occupied molecular orbital (HOMO). Higher 
ionization energies are obtained from deeper molecular 
orbitals. As reference, we show in Fig. the predic- 
tion from LDA (i.e., negative of the Kohn-Sham eigen- 
values for occupied orbitals), although such comparison 
should be made with caution since, from a strict point 
of view, such eigenvalues enter in DFT as mathemati- 
cal entities, with no explicit physical meaning. By ex- 
tending Koopman's theorem to DFT, the HOMO eigen- 
value can be shown to be equal to the first ionization 
potential if the exact functional is used', but such prop- 
erty does not hold for higher ionization potentials". 
Hybrid functional^ and the ones of the exact-exchange 
typ cr^ 2 ^? can give very good predictions for the first ion- 
ization potential. 

In contrast, quasiparticle energies do have physical 
meaning, and Fig. [21 shows definitive consistency be- 
tween theory and experiment. In particular, the first 
ionization energy is predicted to be 9.30 eV within GW/, 
with remarkable agreement with the experimental value 
of 9.3 eV^S^ and also with the prediction from quantum 
chemistry calculations: 9.04 eV— A hybrid B3LYP func- 
tional predict 9.74 eV for this ionization energy— The 
GWo approximation predicts the same quantity to be 
9.88 eV. The difference is due mostly to the vertex term: 
the contribution £/ for the HOMO alone is found to be 
almost 0.8 eV, while LDA exchange-correlation correc- 
tions in the polarizability give a contribution somewhat 
smaller but of opposite sign. For other quasiparticle en- 
ergies, the vertex term is always positive (i.e., it lowers 
the ionization energy compared to the GWo value), and 
no more than 1 eV. 

Ionization energies in the range of 20 to 25 eV are 
found to be systematically underestimated with respect 
to experiment, as shown in Fig. [21 For such deep orbitals, 
we expect to have some loss of numerical accuracy since, 
although the energy dependence of the polarizability n / 
is exactly calculated using Eq. J2J), the summation over 
s is always limited to a finite number of poles. 

The quasiparticle energy for the lowest unoccupied 
molecular orbital (LUMO) is the negative of the elec- 
tron affinity^, and its calculated value is shown in Ta- 
bled Although the anion CgH 6 is unstable, Burrow and 
collaborators^ have been able to perform careful electron 
transmission experiments and identify resonances in the 
spectrum, which where associated to the electron affin- 
ity. Two resonances where found: at energies 1.12 eV 
and 4.82 eV. They match theoretical predictions from the 
GW/ approximation both in energy and spacial distribu- 
tion: e2 U (7r*) and b2 g {^*) respectively. Again, there is a 
discrepancy between GWo an d GW/ quasiparticle ener- 
gies of about 0.5 eV, mostly due to vertex corrections in 
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FIG. 3: Absorption spectrum of benzene, calculated within 
TDLDA (upper panel) and BSE (lower panel). Absorption 
lines were broadened by Gaussian convolution with disper- 
sion 0.15 eV (below 10 eV) and 0.5 eV (above 10 eV). The 
measured spectrum— is shown in dashed lines. In the lower 
panel, the BSE was solved using two approximations for the 
self-energy: GWo (dotted line) and GW/ (solid line). 



the self-energy. 

If measuring electron affinities of benzene is not triv- 
ial, a similar statement can also be made regarding first- 
principles calculations. The anion CgHg is correctly pre- 
dicted within DFT to be unstable, and methods such as 
ASCF (energy variations in self-consistent field) fail to 
predict its electron affinity— Indeed, confining the an- 
ion inside a spherical boundary of radius R and solving 
the Kohn-Sham equations results in an excess energy rel- 
ative to the neutral system that is proportional to Rr 2 . 
This excess energy vanishes in the limit of very large ra- 
dius. The physical picture is that, although the anion is 



TABLE I: Electron affinities in benzene, as predicted within 
GWo, Eq. J22J and GW/, Eq. The outer- 

valence Green's function (OVGF) method— is based on a 
Hartree-Fock expansion of the self-energy, combined with self- 
consistent calculations of self-energy and Green's function. 
All energies in eV. 



LDA 


GWo 


GW/ 


DFT(B3LYP)££ OVGF^ ExppS 


e2u 1-33 


-0.47 


-0.99 


-0.88 -2.605 -1.12 


b 2g -2.45 


-4.49 


-5.05 


-4.82 
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TABLE II: Excitation energies of the lowest-energy neutral excitations in benzene. Energies in eV. 



TDLDA GWp+BSE GW/+BSE TDDFT(B3LYP)^ TDDFT(LHFX)^ CASPT2&2 Exp-^ 



Triplet 
B lu 4.53 


3.59 


3.59 


4.45 


4.27 


3.89 


3.9 


Singlet 














B 2u 5.40 


4.74 


4.86 


5.24 


5.40 


4.84 


5.0 


B la 6.23 


6.08 


6.14 


6.09 


6.12 


6.30 


6.2 


E lu 6.9-7.2 


7.16 


7.23 


6.90 


6.96 


7.03 


6.9 



unstable, it can be detected as a resonance if electrons do 
not remain in their ground state, and the GW approxi- 
mation predicts that resonance as a quasiparticlc orbital. 

As mentioned in the context of Eq. i|32|) . quasiparti- 
cle wave-functions are obtained by numerically diagonal- 
izing the quasiparticlc Hamiltonian, and therefore they 
are different from Kohn-Sham eigenfunctions. For the 
benzene molecule, we observed though that quasiparti- 
clc and Kohn-Sham eigenfunctions overlap by more than 
95 % for all occupied orbitals. This is due to the high 
symmetry of the benzene molecule, which makes most 
off-diagonal matrix elements of the self-energy zero by se- 
lection rules. Similarly the resonant states e-iu and &2g in- 
Table |J have overlaps of 99 % and 63 % respectively. For 
them, overlap is still suppresed by selection rules but the 
existence of many low-energy, unbound orbitals makes 
overlap somewhat more favorable for orbitals well above 
the HOMO-LUMO gap. 

Optical excitations were obtained in two methods: 
both within TDLDA and by solving the BSE. Fig. 
shows the cross section for photoabsorption in benzene. 
It is dominated by a sharp and well pronounced peak at 
around 7 eV, as verified by various calculations^SiiLSS^. 
This peak is the visible component of a tt — tt* complex, 
originated from transitions between the HOMO and the 
quasiparticle state e2 M . The low, flat feature in the 6.0- 
7.0 eV range is due to coupling between tt — tt* transitions 
and vibrational modes^SS, and it is absent in the calcu- 
lated spectra because of the assumed structural rigidity. 
Beyond 10 eV, a number of sharp features in the mea- 
sured spectrum results from transitions involving Ryd- 
berg states£2i§i. The limited numerical accuracy in that 
energy range prevents a detailed identification of such 
transitions in the calculated spectrum. 

Tabic [n] shows a comparison between TDLDA and 
BSE predictions for some excitations in the tt — tt* com- 
plex, together with CASPT2 calculations^ and other 
TDDFT calculationa 52 ! 63 . Although singlet transitions 
B\ u and E\ ul which are the dominant ones in the low- 
energy part of the absorption spectrum, are equally well 
described by both methods, there is a significant blue 
shift of dark transitions B\ u and B\ u within all TDDFT 
approaches presented. With the exception of the E\ u 
transition, excitation energies obtained by solving the 



BSE are typically underestimated with respected to mea- 
sured quantities. The overall deviation between mea- 
sured transition energies and GW predictions is 0.1 to 
0.3 eV, comparable to CASPT2 results^. 

Although excited states can be obtained from either 
TDLDA or BSE by solving suitable eigenvalue problems, 
these two theories provide very different descriptions for 
"electron-hole" correlations. Within TDLDA, such cor- 
relations are included in the exchange-correlation kernel, 
which is static and local in space. Within BSE, they are 
included both in quasiparticle energies (through the self- 
energy operator) and in the kernels K d and K' . And 
these differences manifest themselves in a very intriguing 
way in the tt — tt* complex: both theories agree within 
0.1 eV in the excitation energy of components B\ u and 
Ei u , but they differ by more than 0.5 eV in the excitation 
energy of component £?2 U . 

As observed in Fig. [2] and Tabled the inclusion of ver- 
tex corrections together with a TDLDA polarizability is 
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TABLE III: Electron affinities in naphthalene, Energies in eV. 



LDA GWp GW/ DFT(B3LYP>a QVGF^ Exp-^ 

b lg 1.38 -0.09 -0.60 -0.90 

b 2g 2.19 0.64 0.14 -0.20 -1.30 -0.19 

b Zu 0.76 -1.07 -1.35 -1.67 

a u -1.14 -3.03 -3.44 -3.38 



essential for an accurate prediction of the ionization po- 
tential and electron affinity of benzene. In fact, the GWo 
approximation predicts those quantities to be 9.84 eV 
and -0.54 eV respectively. Defining a "HOMO-LUMO 
gap" as the difference between ionization potential and 
electron affinity, both GWo an d GW/ approaches agree 
in the value of the gap: 10.3 eV. This is somewhat consis- 
tent with the observation by del Sole and collaborators^, 
who conducted a similar analysis in bulk silicon and 
found significant shifts in the valence band maximum 
and conduction band minimum, but only small change 
in the energy gap itself. Due to the cancellation of ver- 
tex contributions, the energy position of the 7r — 7r* line 
is not significantly affected by the choice of self-energy 
approximation in the solution of the BSE, as shown in 
the lower panel of Fig. [3J 

Naphthalene, CioHg is the second smallest oligoacene. 
For this molecule, we used the following parameters in 
the solution of the Kohn-Sham equation: grid spacing 
0.4 a.u., boundary radius 20 a.u., and experimental val- 
ues of bond length. Ionization potentials associated to all 
occupied 7r orbitals and some a orbitals were reported by 
Dewar and Worlej^, with the first one being 8.11 eV— . 
Some of the theoretical calculations of the first ionization 
potential arc: 7.85 cV (OVGF method^); and 7.59 eV 
(hybrid B3LYP/DFT 68 ). Predictions from the GW/ and 
GWo approximations are 8.15 cV and 8.69 eV respec- 
tively. Fig. 01 presents a comparison between predictions 
from the GW/ and GWo approximations for all other po- 
tentials. Although both are acceptable, the GW/ approx- 
imation is clearly superior to GWo, giving a maximum 
deviation from experiment of no more than 0.3 eV. Ver- 
tex corrections are of the order of 0.6 cV. We should point 
out that the magnitude of self-energy corrections for ir or- 
bitals (experimental ionization energy ranging from 8 to 
13 eV— ) is around 2.4 eV, whereas the self-energy cor- 
rections for a orbitals (experimental ionization energy 
ranging from 13 to 19 eV— ) is bigger: around 3.6 eV. 
Since a orbitals are typically deeper in energy than 7r, 
they are expected to have self-energy corrections larger 
in magnitude. 

There has been some debate in the literature concern- 
ing the stability of the anion CioHg— . Early electron 
capture experiments have observed a stable anion with 
electron affinity in the range of 0.15 eV— . In contrast, 
electron transmission spectroscopy measurement have in- 
dicated the ion to be unstable, with negative electron 



affinity of -0.2 eV— . The anion is predicted by GW/ 
to be stable, with electron affinity of 0.14 eV. Table ITTT1 
presents a comparison between the calculated electron 
affinities and the experimental data^. The fact that the 
anion has such small binding energy makes any stabil- 
ity analysis, cither from theory or experiment, very dif- 
ficult. As observed in benzene, there is a number of 
low-energy orbitals around the LUMO, most of them 
highly delocalized. But some of those orbitals are lo- 
calized around the molecule, and they are responsible 
for the resonant states detected in electron transmission 
spectroscopy measurements^. For the resonant states, 
we find fair agreement in terms of orbital character and 
corresponding electron affinity, although there is a dis- 
crepancy of typically 0.2 to 0.3 eV in the last quantity. 

The absorption spectrum of naphthalene is found to 
have a sharp and well-pronounced line at 6.0 eV. This 
line is originated from a B\ u transition which is optically 
active for light polarization along the line that joins the 
centers of the two aromatic rings in the molecule, called 
"long molecular axis" . A second, lower absorption line is 
located at 4.3 eV and it has B 2u character, being highly 
active for polarization on the plane of the molecule but 
perpendicular to the long molecular axis. TDLDA and 
BSE give very similar values for the energy position of 
both lines, although they differ by 0.5 eV in the predicted 
value of the first excited state, a dark B\ u singlet excita- 
tion. A comparison between the theoretical calculations 
and experimental data is shown in Tabic llVl 

In summary, the absorption spectra of benzene and 
naphthalene are dominated by a sharp 7r — 7r* transi- 
tion, found in both TDLDA and GW+BSE methodolo- 
gies. Considering that this transition is composed of hy- 
bridized carbon-p orbitals, it is expected to be very lo- 
calized in space. We do not discard the possibility of 
strong spacial localization being a major reason for such 
good agreement between TDLDA and GW+BSE, despite 
the conceptual simplicity of the former. The good agree- 
ment between predicted and measured ionization poten- 
tials of benzene and naphthalene is equally remarkable. 
There has been a very limited number of cases where such 
quantities were calculated within the framework of the 
GW approximationii2i2iiLi2iiS, and the present results 
provide the first benchmarks for ionization potentials in 
oligoacenes. 
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TABLE IV: Excitation energies of the lowest-energy neutral excitations in naphthalene. Spin-singlet states only. Energies in 
eV. 



TDLDA GWo+BSE GW/+BSE TDDFT(B3LYP)^ CASPT2^ Exp^ 





4.40 


3.93 


4.02 


4.21 


4.03 


3.97 


B211 


4.32 


4.28 


4.34 


4.12 


4.56 


4.45 


Bin 


5.84 


6.04 


6.12 


5.69 


5.54 


5.89 


B2u 


6.13 


6.08 


6.16 


5.90 


5.93 


6.14 



B. Silicon clusters 

Small semiconductor clusters containing a few tens of 
atoms and with various chemical compositions have been 
mass produced and characterized in terms of their op- 
tical properties and structural configuration, although 
the synthesis of small silicon clusters remains a chal- 
lenging tasl*i^2i*£^. From a theoretical point of view, 
optical properties of small silicon clusters are interest- 
ing for one aspect: the linear optical spectrum of bulk 
silicon is known to be extremely well predicted by a 
first-principles approach based on the BSFj£*£, whereas 
TDLDA fails to predict both the optical gap and exci- 
tonic effects^iiSSi. On the other hand, optical excita- 
tions of S1H4 and Si2H 6 , the two smallest "clusters", are 
correctly described within TDLDA2&. Based on these two 
facts, it is natural to investigate what limits the validity 
of each theory and where the crossover size is, beyond 
which one theory becomes superior to the other. Such an 
analysis is challenging because of two major factors: the 
limited number of useful experiments, and the complexity 
of some numerical implementations of the BSE method. 
Recently, a study based on model dielectric screening in 
silicon clusters have addressed the issue, with inconclu- 
sive results^. Since the present implementation of the 
BSE method is particularly efficient for finite systems, 
we are able to compare both TDLDA and BSE meth- 
ods in a wide range of cluster sizes and with a minimum 
of ad hoc assumptions. As we discuss below, there is un- 
fortunately one additional difficulty in such comparisons: 
there is no direct relationship between the optical spec- 
trum of a nanostructure (quantified by the absorption 
cross section, for instance) and the optical spectrum of a 
macroscopic solid (quantified by the macroscopic dielec- 
tric function). 

In general, the surface of synthesized silicon clusters 
is covered by passivating particles and, for the smallest 
ones, they may undergo significant reconstruction^. Our 
approach is to construct clusters from fragments of crys- 
talline silicon and adjust the number of atoms so that 
the fragment has surface as spherical as possible. Dan- 
gling bonds on the surface are passivated by attaching 
hydrogen atoms. We do not consider surface reconstruc- 
tion. Although the class of tetrahedral clusters we discuss 
here may not be the most stable ones in the size range 
of to 2 nm, they are expected to be the most stable 



ones when the cluster size is very large. Here, we ana- 
lyze the clusters SiH 4 , Si 5 Hi 2 , Sii Hi 6 , Sii4H 2 o, Si 2 9H 36 , 
S135H36, Si47H 60 , Si7iH 8 4, SiggHioo, and Sii47Hioo- Since 
the presence of hydrogen atoms typically requires dense 
grids, we use grid spacings ranging from 0.5 a.u. (S1H4) 
to 0.7 a.u. (Sii47Hioo) in the solution of the Kohn-Sham 
equations. 

For each cluster, the quasiparticle Hamiltonian is re- 
orthogonalizcd in the Kohn-Sham basis. As emphasized 
in the literature&ii, the self-energy operator is not diago- 
nal, and Kohn-Sham cigenfunctions are not good approx- 
imations for quasiparticle wave- functions, particularly for 
the low-energy unoccupied orbitals in the smallest clus- 
ters. There are two main reasons for this phenomenon: 
LDA wrongly predicts some of these unoccupied states 
to be bound, whereas quasiparticle orbitals are often un- 
bound and more extended; and the large density of or- 
bitals within a few eV from the LUMO makes the oc- 
currence of non-negligible off-diagonal matrix elements of 
the operator in Eq. I).'.i2[l very frequent. The Si35H36 clus- 
ter is a typical example. There, overlap between Kohn- 
Sham and quasi-particle wave- functions is observed to 
be greater than 95% for most occupied states and for the 
first few unoccupied ones, but it reduces to values rang- 
ing from 30% to 90% for the higher-energy unoccupied 
states. 

One important aspect involving self-energy corrections 
is whether they can be modeled by a "scissors" operator 
or not. Frequently, the quasiparticle band structure of 
solids differs from the one predicted within LDA by a 
rigid upward shift of conduction bands with respect to 
valence bands2<4, which justifies the scissors approxima- 
tion. Scissors operators with energy dependence can also 
improve band widths^. But we find this procedure to 
be inaccurate even for moderately large clusters. Fig. 
shows diagonal matrix elements of the self-energy opera- 
tor as function of the LDA energy eigenvalue for all oc- 
cupied orbitals and some unoccupied orbitals of Si35H 3 6. 
For deep occupied orbitals, these matrix elements follow 
a quasi-continuous distribution, with smooth curvature, 
so a scissors operator could be defined. But the matrix 
elements for orbitals in the vicinity of the energy gap do 
not have a well-defined dependence with respect to LDA 
eigenvalues. On the contrary, they have strong orbital 
dependence. This fact, together with the existence of 
large off-diagonal matrix elements, makes any attempt 
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at constructing a scissors operator very difficult for these 
clusters. 

Fig. [3] shows the dependence of electron affinity and 
first ionization potential with respect to the cluster diam- 
eter. To our knowledge, the ionization potential has not 
been measured for any of the studied clusters except for 
S1H4, for which the experimental value is 12.6 eV—. For 
this system, the prediction from GW/ is 12.5 eV. Starting 
from SiH4, the first ionization potential decreases mono- 
tonically as the cluster size increases. On the other hand, 
the electron affinity remains very small, with little depen- 
dence with respect to cluster size. As a result, the elec- 
tronic gap (defined as the difference between ionization 
potential and electron affinity) decreases continuously for 
larger and larger clusters. This is a size effect which has 
been observed before for small clusters^. In order to un- 
derstand the physics behind it, we model the cluster as an 
electron sphere with homogeneous density and radius R, 
keeping the density constant and proportional to N/R 3 . 
The energy of the HOMO is found by integrating the 
density of states up to the number of electrons N, result- 
ing m E F = const, x N 2 / 3 R" 2 . The energy difference 
between this orbital and the next one is approximately 
AE = const, x N~ 1 / 3 R~ 2 = const, x i?~ 3 , which is in- 
versely proportional to the volume of the cluster. A more 
elaborate analysis predicts a power law R~ 2 instead of 
R~ 
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FIG. 5: First ionization potential and electron affinity of pas- 
sivated silicon clusters, calculated within the GW/ approxi- 
mation (solid lines) and ASCF (dotted lines). Experimental 
data from Refi— . ASCF results include spin-polarization ef- 
fects. 



interaction effect, present in ASCF calculations, and the 
incorrect asymptotic behavior of the LDA functionali* 7 ^. 

Fig. [7| shows the absorption cross section for clusters 
SifLi to Si47H6o- The cross section also has a characteris- 
tic behavior as function of cluster size: for small clusters, 
it has a small number of well-defined peaks but, as size 
increases, their intensity decreases in the low end of the 
energy axis. For large clusters, the absorption cross sec- 
tion shows a smooth and featureless profile, with onset 
at a low energy value and continuous increase towards 
higher values of energy. Also, TDLDA and BSE spectra 
are significantly different for small clusters. In particular, 
TDLDA predicts as-tp line at around 8.2 eV, whereas 
the s — > p line in BSE is positioned at 9.4 eV. Previous 
BSE calculations have predicted 9.0 eV-£ and 9.16 e\£ 
for the same quantity. The experimental value has been 
reported to be around 8.8 eV— . For bigger clusters, the 
differences between TDLDA and BSE tend to disappear. 

As pointed out recently^, the behavior of the absorp- 
tion cross section for large clusters is dictated by the ex- 
change kernel K£ cv , c ,, which is present in both TDLDA 
and BSE formalisms. This term introduces a long-range 
interaction that increases in strength as function of clus- 
ter size and ultimately dominates both the exchange- 
correlation kernel in the TDLDA eigenvalue equation and 
the direct+vertex kernel in the BSE. As a result, oscil- 
lator strength is redistributed among the excited states 
of the system and the absorption spectrum peaks at fre- 
quencies well above the optical range. Differences in the 
detailed treatment of electron-hole correlations, intrinsic 
to both TDLDA and BSE, are expected to be less and 
less important for increasingly large clusters. 

As the cluster size increases, more transitions become 
allowed and the absorption cross section is expected to 
evolve into a well-pronounced and broad peak at the 
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We observe a systematic discrepancy in the first ion- 
ization potential as predicted within GW/ and ASCF, 
the difference being 0.6 to 0.9 eV throughout the range 
of studied cluster sizes. In contrast, the electron affinity 
shows better agreement between these two theories. Pos- 
sible explanations for this behavior are the spurious self- 



FIG. 6: Diagonal matrix elements of the operator E — V xc , 
c.f. Eq. 132H . evaluated in the basis of Kohn-Sham eigenfunc- 
tions. Matrix elements for occupied (unoccupied) orbitals are 
represented by crosses ("plus" signs). Self-energy is evalu- 
ated within the GW / approximation and at the extrapolated 
quasiparticle energy (see text). 
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tides in very low concentration, the imaginary part of 
the dielectric function is proportional to the absorption 
cross section of each particle according to the expression 
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FIG. 7: Cross section for photoabsorption in passivated sil- 
icon clusters, calculated within TDLDA (dashed lines) and 
BSE (solid lines). The self-energy operator used in BSE is 
calculated within the GW/ approximation. The cross section 
is normalized by the number of silicon atoms in each cluster. 
Absorption lines were broadened by a Gaussian convolution 
with dispersion 0.1 eV. 



plasmon-pole frequency. Indeed, the formalism in Sec- 
tions and ^TCl can be applied directly to solids. In 
solids, Eq. (|llf) reduces to a function proportional to the 
imaginary part of the inverse dielectric function, mea- 
surable by electron energy loss spectroscopy, EELS. 01- 
evano and ReiningZi have shown that the random-phase 
approximation correctly predicts the essential features of 
the EELS of bulk silicon, and including corrections at 
BSE level does not lead to substantial improvements. As 
expected, the spectra of the largest clusters studied, de- 
picted in the lower panels of Fig. show a smooth profile 
typical of the low-energy end of the plasmon peak. 

The blue shift can be discussed in terms of classical 
electrodynamics. Wc consider a material composed of 
several identical, spherical particles in suspension inside 
an optically inert background. According to the Mie the- 
ory (or effective medium theory)^2&, the (macroscopic) 
dielectric function of the medium is size-dependent, and 
there are three important regimes: particles much smaller 
than the wavelength of light; particles of size compara- 
ble to the wavelength of light; and particles much bigger 
than the wavelength of light. In the regime of small par- 



£2 = nkuJI'K 



(41) 



where A is the wavelength of light and n is the concentra- 
tion. The two functions then share the same energy de- 
pendence. This is a special case of the Clausius-Mossotti 
relation when n — ► 0. In denser medium, the proportion- 
ality above no longer holds and the energy dependence 
of £2 is now given by a non-linear Clausius-Mossotti rela- 
tion. In the extreme situation of highly packed particles, 
with no interstices, the function £2 should resemble the 
bulk dielectric function. As the particles condense, any 
similarity between the energy dependencies of £2 and a 
is lost. The qualitative aspects of this analysis still hold 
for more complicated cases, such as non-spherical and/or 
inhomogeneous particles. Some predictions of the Mie 
theory for the optical properties of nanostructures have 
been discussed in recent first-principles calculations: the 
emergence of a Mie plasmon in nanoclusters^; and the 
dependence of dielectric function with respect to concen- 
tration in periodic arrangements of nanowires^. 

In Fig. [7| the parameter is particle size rather than 
concentration. From the calculated spectra, one can re- 
cover the bulk limit by increasing the concentration or 
the size of particles. In the first case, the Mie theory pre- 
dicts a red shift of £2 with respect to a. The second case 
should certainly involve a similar phenomenon but it re- 
quires a more elaborate analysis, because of the crossover 
between particle size and wavelength. The crossover does 
not necessarily affect the absorption cross section because 
light is not an essential element in absorption measure- 
ments (for instance, one can use electron beams instead 
of photons as probe). But it should affect the dielectric 
function because measurements of dielectric function nec- 
essarily involve interaction with photons. 

Although the behavior of the absorption cross section 
is well understood, we conclude that spectra of large clus- 
ters are not directly related to the macroscopic dielectric 
function. Whereas excitonic effects and optical band gaps 
are easily recognized in the macroscopic dielectric func- 
tion, the same docs not hold for the inverse dielectric 
function. In principle, these two functions are linked by 
simple inversioniL^i, but the delicate features of both 
functions can be lost if numerical inversion is attempted 
without careful control of numerical accuracy. For solids, 
numerical inversion is avoided by replacing the BSE for 
the polarizability operator with a similar equation that, 
when solved, provides the macroscopic dielectric function 
directlyiSi. For confined systems, no such prescription 
is known, and therefore the question of how to compute 
an absorption spectrum that, in the bulk limit, reduces 
unambiguously to the macroscopic dielectric function re- 
mains unanswered. 

The problem of an indeterminate crossover and ill- 
defined bulk limit in the absorption spectrum has two 
important consequences. The first one is that the iden- 
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tification of excitonic effects in large silicon clusters be- 
comes a non-trivial task, even if the cluster is sufficiently 
large so that the very concept of "exciton" applies. The 
second is that an optical gap cannot be easily extracted 
from the absorption cross section of large clusters. As 
seen in Fig. the onset of absorption for the large clus- 
ters is not well defined. Often, low-energy excitations in 
the cluster have extremely small oscillator strength and 
they should not be used to define a gap. Presently, one 
tentative definition of the optical gap is via a threshold 
approach: the optical gap Eq is such that the integrated 
cross section up to Eq is a fixed fraction of the total 
integrated cross section^: 



pEa poo 

/ a(E)dE = p a(E)dE , 
Jo Jo 



(42) 



where p is a small fraction. This prescription is motivated 
by experiment: often, the onset of absorption cross sec- 
tion cannot be distinguished due to limitations in exper- 
imental resolution, and a practical definition of optical 
gap becomes important. An obvious drawback of Eq. 
(|42[) is that the gap may be sensitive to the choice of 
p. In practice, this prescription has been observed to be 
robust. 

In Fig. |H| we show the minimum gap and optical gap 
as predicted within TDLDA and BSE. The minimum gap 
is defined as the energy of the lowest excited state, irre- 
spective of oscillator strength, whereas the optical gap is 
defined from Eq. P^l with p = 2 x 10~ 4 . As expected, 
both gaps have strong dependence with respect to cluster 
size. In addition, the BSE gaps are typically larger than 
the corresponding TDLDA ones, the difference slowly 
decreasing for larger clusters and almost vanishing for 
Sii47Hioo- Fig. [8] also shows that the difference between 
TDLDA minimum and optical gaps has a slow tendency 
towards increasing with cluster size, but the same is not 
observed in the BSE gaps. The reason for this behavior 
is that low-energy excitations obtained within BSE have 
larger oscillator strengths than the corresponding ones in 
TDLDA. Discrepancies between the TDLDA optical gap 
in Fig. [5] and previous work2£ is due to different choices 
for the value of p. 

The dependence of the minimum gap with respect to 
cluster size has also been discussed in the context of 
quantum Monte-Carlo (QMC) calculations^. QMC and 
GW+BSE share a common link: both of them are many- 
body theories, and therefore they both include correla- 
tion effects not present in TDLDA. A signature of these 
effects can be recognized in the calculated minimum gap: 
QMC and GW+BSE give similar values for that quan- 
tity, both of them overestimated with respect to TDLDA 
(c.f. Fig. [S]and RefiS). For large clusters, there is a ten- 
dency for the QMC gap to fall above the GW+BSE gap. 
For instance, the GW+BSE minimum gap for Sii47Hioo 
is found to be 3.3 cV, whereas QMC predicts a gap of 
around 4 eV— . The corresponding gap within TDLDA 
is found to be 2.5 eV. 
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FIG. 8: Minimum gap (dashed lines and open symbols) and 
optical gap (solid lines and filled symbols) for passivated sil- 
icon clusters. Circles denote TDLDA gaps. Triangles denote 
GW/+BSE gaps. The optical gap is obtained by integrating 
the oscillator strength up to p = 1 x 10~ 4 . The minimum gap 
is defined as the lowest spin-singlet excitation energy (p — 0) 
in the system. Experimental data from Refs^i*— . 



C. F-center defects in LiCl 

The formalism presented in Section lll Bl can also be ap- 
plied to extended systems, once the appropriate bound- 
ary conditions are used in the solution of the Kohn-Sham 
equations^. In fact, the boundary conditions imposed on 
electron wave-functions determine the ones for the po- 
larizability operator II, through Eqs. J2J), ©, iffil . and 
(|4Ufl . In periodic systems, optical absorption is quanti- 
fied by the absorption coefficient, or by the imaginary 
part of the macroscopic dielectric function. Here, we cal- 
culate this dielectric function following an approach due 
to Hank o 7 ! 81 . The basis of this approach is in defining 
a truncated Coulomb potential V in the construction of 
the exchange kernel K x : 



K x — > K x 

vcv'c' vcv'c' 



Aire 2 (v\X ■ v|c) (c'|A • v\v') 

AfcV^ e // £ c £ v S c ' S v ' 



where A • v is the velocity operator projected along some 
direction A. The volume of the periodic cell, V ce u and 
number of k-points, Nk, are included as normalization 
factors. The kernel K x in the above equation does not 
have direct physical meaning, being rather an auxiliary 
quantity, which differs from the original K x by the ab- 
sence of a long-range Coulomb interaction. By using this 
truncated kernel, Hankc has shown^ithat a polarizability 
operator can be obtained and, from that, the macroscopic 
dielectric function according to 
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(43) 

where p s and uj s can be found by diagonalizing cither 
TDLDA-like or BSE-like eigenvalue equationo 7 i 73 i 81 . In 
addition, lattice periodicity allows us to use a wave- vector 
q as additional quantum number in the polarizability IT. 
Each choice of q leads to a different eigenvalue equation 
and the resulting partial polarizabilities are summed up. 

Here, we are interested in analyzing the optical spec- 
trum of a F center embedded in an otherwise ordered LiCl 
crystal. F centers are halogen vacancies found in alkali- 
halide crystals, and they are named after the bright, vis- 
ible coloration they induce. Alkali-halide crystals are 
well-known for their strong ionicity and wide energy gap, 
which makes them transparent to visible light. Halogen 
vacancies induce localized electronic states inside the gap, 
and electronic transitions involving such localized states 
are the source of the bright color in F-center-rich crystals. 
The nature of F centers has been studied experimentally 
using various techniques^. 

From a theoretical point of view, F centers are chal- 
lenging because the exact location of defect states in- 
side the gap is not easily obtained without experimen- 
tal input. Effective-mass models are not suitable due 
to the extremely large band gap, high effective masses, 
and small dielectric constant. Also, DFT suffers draw- 
backs due to the intrinsically underestimated band gap. 
Within the GW approximation, not only the band gap is 
correctly predicted, but defect states around the F center 
in LiCl can also be located^. Here, we go beyond and 
determine optical excitations and the dielectric function 
within BSE, and compare the results with the TDLDA 
prediction and with experiments. 

Defect levels in LiCl are extremely localized, which 
simplifies considerably the numerical work. We simu- 
late the defect by using a cubic supercell with 32 lithium 
and 32 chlorine atoms. The lattice parameter is taken 
from experiment. Starting from the ordered structure, 
one chlorine atom is removed and the Kohn-Sham equa- 
tions arc solved in real space with grid spacing 0.3 a.u. 
All atoms inside the periodic cell are allowed to move 
so as to minimize the total energy. Only atoms in the 
first lithium shell showed significant relaxation, with dis- 
placement of 0.05 a.u. and relaxation energy 30 meV. As 
convergence test, the procedure was repeated in a much 
bigger cell, containing 216 atoms of each type (excluding 
the vacancy). Movement of atoms in the first shell was 
around 0.02 a.u. and relaxation energy less than 10 meV, 
thus showing the degree of localization of the defect. In 
the subsequent work, we used the 32Li+32Cl+vacancy 
cell. Since the presence of a vacancy breaks translational 
invariance, we reduced the Brillouin zone to the T point 
only, thus removing complications with Bloch functions 
and integrations over the Brillouin zone£4. 

Quasiparticle energies for the band edges and two de- 



TABLE V: Electron band structure of the F center in LiCl, 
with defect levels. All calculated energies are given with re- 
spect to the valence band maximum as predicted by GW/, in 
eV. Results obtained from a 2x2x2 cubic supercell, containing 
63 atoms. Defect level 2p has a spin splitting of 0.4 eV. Re- 
sults in column "GW bulk" were obtained by Hybertsen and 
Louie- for bulk LiCl. 
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feet levels are shown in Tabled The choice of supercell 
allows us to recognize easily the band edges at points 
r, L, and X by observing the symmetry properties of 
the calculated single-particle orbitals. Self-energy cor- 
rections in the neutral F center require the explicit in- 
clusion of spin degrees of freedom: the Is state is par- 
tially occupied with only one electron. Therefore, the 
self-energy operator was evaluated separately for the two 
possible spin configurations and the summation over oc- 
cupied states in Eq. I|23|) performed over all occupied 
states for each spin configuration. A similar procedure 
was used for the evaluation of the correlation and vertex 
parts of the self-energy. Although the self-energy oper- 
ator is spin-dependent, most quasiparticle levels remain 
spin-degenerate, which is expected since there is no in- 
trinsic source of spin polarization in the system. Table 
IVl shows that, within DFT-LDA, the Is state is posi- 
tioned inside the gap, but the 2p triplet state is located 
within the conduction band. Inclusion of self-energy cor- 
rections keeps the 2p triplet above the conduction band 
maximum (CBM). This phenomenon has been observed 
before^. Although mixed with extended states, the 2p 
triplet remains fairly localized within the vacancy. In- 
deed, 45 % of its probability distribution is located within 
a distance of Lq or less from the vacancy, where L$ is the 
Li-Cl interatomic distance. For comparison, 74 % of the 
Is probability distribution is located within the same re- 
gion. 

Optical excitations in the system are calculated in both 
TDLDA and BSE frameworks. In both cases, spin po- 
larization is included explicitly in the construction of the 
electron-hole kernel. Although not essential at TDLDA 
level, the inclusion of spin polarization is important in 
the BSE because of the existence of the Is partially oc- 
cupied orbital. Fig. [!|] shows the imaginary part of the 
dielectric function for the cubic supercell containing 63 
atoms. Although this crystal has an electronic band gap 
of 9.4 eV, the BSE spectrum is dominated by a wide dou- 
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FIG. 9: Imaginary part of the macroscopic dielectric function, 
£2 for a F center embedded in LiCl. Solid and dashed lines cor- 
respond to GW/+BSE and TDLDA calculations respectively. 
The self-energy operator used in BSE is calculated within the 
GW/ approximation. The position of the Is — > 2p transi- 
tion line and electronic bands gaps within LDA and GW / are 
shown with arrows. A Gaussian convolution with dispersion 
0.2 eV was used to smear out the absorption lines. 



ble peak just below 9 eV. Being located below the band 
gap, this peak has the signature of excitonic effects^. The 
limited supercell size compromises resolution of the exci- 
tonic double peak and causes a redshift of 0.2 to 0.4 eV, 
although the overall features are consistent with calcu- 
lations for a clean LiCl crystal. Below the excitonic 
peak, we can see one absorption line corresponding to the 
Is — > 2p transition. In the real material, the strength of 
this line depends on temperature and density of defects, 
and so the calculated strength is somewhat arbitrary. 

The absorption spectrum obtained within TDLDA is 
characteristic for the absence of the exciton peak. In 
fact, it has an onset at the DFT-LDA band gap, around 
6.3 eV, followed by a featureless rise towards higher en- 
ergies. The three wide peaks in Fig. located at 6.4, 
8.3 and 9.7 eV are due to extremely limited resolution in 
the supercell. By including more atoms in the supercell, 
these peaks are expected to merge into a smooth func- 
tion. The lack of excitonic effects in TDLDA has been 
reported in the Iiteraturo£ii2i2i. 

Surprisingly, TDLDA does predict a very accurate 
value for the Is — > 2p transition energy: 3.10 eV, com- 
pared to 3.04 eV (BSE) and the measured value of 
3.25 eV to 3.3 e V^'?7i9?,??. We believe this is related 
to the strong localization of both Is and 2p orbitals. As 
discussed before, these orbitals are confined to within 1 
or 2 atom sites from the vacancy. To some extent, the 
defect behaves as a confined electronic system, separated 
from the crystalline background. For transitions internal 
to this defect, such as the Is — > 2p, the LDA exchange- 
correlation kernel correctly describes electron-hole inter- 
actions. But transitions involving the conduction and/or 
valence bands are not so well described because the LDA 
kernel lacks the non-locality present in electron-hole in- 



teractions involving extended orbitals. The BSE spec- 
trum also shows a low feature at around 5.8 eV, which 
can be assigned to the Is — > L lineS^. 

IV. CONCLUSION 

We discussed an implementation of Green's functions 
methods developed in the space of single-particle transi- 
tions. Contrary to alternative approaches, this procedure 
does not make use of Fourier analysis, and thus can be 
applied directly to confined systems, where it is partic- 
ularly efficient. As an example, we are able to do a full 
GW+BSE calculation of the benzene molecule using a 
1.7 GHz IBM Power4 machine in single- process mode in 
280 minutes. This task requires no more than 600 MB of 
CPU memory and the only input required is the geom- 
etry of the molecule. Taking advantage of the extreme 
numerical efficiency of the current implementation, we 
are able to perform ab initio GW+BSE calculations of 
silicon clusters containing more than one hundred atoms. 

Electronic screening is included within the TDLDA, 
and a corresponding vertex is included for consistency in 
the diagrammatic expansion. The added terms (TDLDA 
vertex + TDLDA screening) are found to be essential 
for accurate predictions of electron affinity and ioniza- 
tion potentials of benzene and naphthalene. This is in 
contrast with previous GW calculations which assume 
screening at RPA level only but nevertheless provide ac- 
curate ionization potentials for small molecules&2iiLiS. 
The explanation for this apparent contradiction may be 
that the popularly used plasmon-pole models carry infor- 
mation about the exact many-body screening and there- 
fore corrects deficiencies of the RPA. 

Besides oligoacencs, the current implementation is 
used to investigate the absorption cross section of sili- 
con clusters by solving the BSE. We conclude that, while 
the TDLDA and BSE differ markedly for small clusters, 
they agree in the broad features of the spectrum for large 
clusters. In particular, the dependence of excitation en- 
ergy as function of cluster size is similar for both theories. 
The present results are consistent with QMC calculations 
for the minimum gap of silicon clusters. Excitation en- 
ergies, ionization potentials and electron affinities calcu- 
lated within the present method are consistent with ex- 
perimental data to within a fraction of eV, comparable 
to chemical accuracy^!. In LiCl, we show how the exis- 
tence of a F center affects the energy-resolved dielectric 
function. We also show that TDLDA predicts correctly 
the position of the Is — > 2p, despite producing a poor 
description of excitonic effects. 
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APPENDIX A: EVALUATION OF ENERGY 
INTEGRAL IN EQ. $22$ 

The derivation of Eqs. P, (23} and from Eq. 
(|22[1 follows from standard integration over poles^l. For 

completeness, we present here the major steps. We start 
by defining separate exchange and correlation contribu- 
tions in Eq. 



mmw) = m x \j') + ois c (W> (ai) 



with 



xGiv^E' -E)i Pj ,(v 2 )V(v 1 ,v 2 ) , (A2) 



APPENDIX B: STATIC REMAINDER IN 
SELF-ENERGY 

In most cases, the summation over n in Eqs. (|24|l and 
(131) I has slow convergence. This behavior is not unique 
to the present implementation and it has also been ob- 
served in self-energy calculations when the dielectric ma- 
trix is explicitly computed in reciprocal spaca2*2*2£. In 
that case, converged self-energies were found to require 
in excess of 100 unoccupied bands in the single-particle 
Green's function. The convergence rate can be accel- 
erated by truncating this summation at some point and 
evaluating the remainder within the COHSEX (Coulomb 
hole + screened exchange) approximation^. The COH- 
SEX approximation is essentially a static limit of Eq. 
(I22|l . Assuming that screening is instantaneous, the po- 
larizability IIo and the potential Wq become proportional 
to a (5-function in time and constant in energy. The in- 
tegral over energy in Eq. (|22|l can then be replaced by 
a sum over poles of the single-particle Green's function. 
In our implementation, one can recover the COHSEX 
approximation by imposing u> s >> \E — e n \ in Eq. i|24|> : 



and 

(j\X c (E>)\f) = fdr 1 dr 2i p j (r 1 )if% e~ iE0+ 

xG( ri ,r 2 ;£'-£)^(r 2 ) x 
[/dr3dr4V(r 1 ,r3)no(r3 ! r 4 ; J B)y(r4 ! r 2 )] . (A3) 

We assume a quasiparticle approximation for the one- 
electron Green's function, 



U\2c(E)\3')\ 



COHSEX 



(Bl) 



The last summation on n is done over all Kohn-Sham 
cigenstates. We evaluate it exactly by using a complete- 
ness relation. Using Eqs. ©, © and lj2"5)l. we ob- 
tain: 



G(v u r 2 ;E) 



^ E - e n + ir] n W 



(A4) 



In the exchange term, Eq. I|A2|I . the only poles present 
are the ones originated from the Green's function. There- 
fore, the energy integration is easily replaced with a sum- 
mation over occupied states only, resulting in Eq. (|23|l . 
For the correlation term, Eq. i|A3|) , we assume a polariz- 
ability operator given by Eq. Now, both G and Ho 
contribute with poles below the real energy axis. Col- 
lecting them, one arrives at Eq. I|24|) . 



0'|£c(W)l 



COHSEX 



where 



= 4V occ V v ™J V nj' 
i/dr^-(r)lF(r)^(r) (B2) 



W(r) = - I dr' / dr"U(r,r')n (r',r";i; = 0)U(r",r) . 



Eqs. (|B1|I and i|B2|) can be used to determine the effect 
of truncating the sum over n at a value N >> n occ . The 
remainder then is 



RC iA N ) = \ J dr^(r)^(r) J dr' J dr"U(r, r')n (r', r"; 0)U(r", r) + 2 £ ]T 



V S V S ; 

nj nf 
UJ S 



(B3) 



r 



Although the COHSEX approximation has a level of with 
accuracy lower than the full GW method, the conver- 
gence behavior is similar in both approximations if N is 
chosen sufficiently high. Eq. I|24|) can then be replaced 
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N 



V-V; 

nj Y nj' 



(B4) 



N 



Along similar lines, the vertex term, Eq. H30|) is rewrit- 
ten as with 



(B5) 



J 



rUn) 



i / tow (r)<Pj> (r) / dr [V(r, r')IL f (r>, r; 0)/(r) 

I 



f(r)U f (r,r';0)V(v',T)} 



(B6) 



As an example, we have computed self-energy correc- 
tions in benzene including and not including the static 
remainder. With a static remainder, the first ionization 
energy decreases by 20 meV when N increases from 256 
to 512, resulting in ionization potential of 9.30 eV for the 
larger value of N. Without static remainder, this ion- 
ization energy increases by 133 meV between the same 
choices of N, and its value for N = 512 is 8.27 eV, which 
is still far from convergence by as much as 1 eV. 

In molecules and clusters, high-energy virtual states 
are expected to be very delocalizcd, and therefore sensi- 
tive to the choice of boundary conditions. In spite of that, 
the use of confined-system boundary conditions (where 



all wave-functions are required to vanish outside some 
spherical enclosure) is still justified. The reason is that, 
as shown in Eq. (|24|l . only the overlap between high- 
energy and low-energy states is relevant for self-energy 
calculations. The detailed features of high-energy states 
in the vacuum region, away from the atoms, are not very 
important. In addition, the contribution of virtual states 
in the summations over n decreases as one goes to higher 
and higher states. Nevertheless, the size of the confining 
region should always be tested against convergence, so 
that the shape of virtual states in the vicinity of atoms 
is correctly described. 
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As written, Eq. (J5| corresponds to a time-ordered func- 
tion, but it differs from the causal (measurable) response 
function only for negative values of energy. 
Confirming earlier work-, we have observed that the 
Tamm-Dancoff approximation (i.e., neglecting mixing be- 
tween absorption and emission contributions) is accept- 
able when calculating excitation energies in the GW+BSE 
framework. This contrasts with the scenario in TDDFT, 
where mixing is usually importanlA— ■ , due to the spe- 
cific properties of the TDDFT kernels. The impact of the 
Tamm-Dancoff approximation in the energy loss spectrum 
of silicon has been analyzed by Olevano and Reining—. 
The value of 7.4 eV reported in Refi— ■ for the peak position 
has been recently revised to 7.1 eV in Ref— ■. 
The inaccurate optical gap reflects the gap underestima- 
tion inherent to DFT, and it can be corrected by an ad hoc 
scissors operator. The lack of excitonic effects should be 
attributed to the specific functional employed in TDLDA, 
since other functionals were shown to correctly predict the 
enhancement of the E\ peak relative to i?2— 



